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Abstract 

We show that the number Z of g-edge-colourings of a simple regular graph of degree q is de- 
ducible from functions describing dimers on the same graph, viz. the dimer generating function or 
equivalently the set of connected dimer correlation functions. Using this relationship to the dimer 
problem, we derive fermionic representations for Z in terms of Grassmann integrals with quartic ac- 
tions. Expressions are given for planar graphs and for nonplanar graphs embeddable (without edge 
crossings) on a torus. We discuss exact numerical evaluations of the Grassmann integrals using an 
algorithm by Creutz, and present an application to the 4-edge-colouring problem on toroidal square 
lattices, comparing the results to numerical transfer matrix calculations and a previous Bethe ansatz 
study. We also show that for the square, honeycomb, 3-12, and one-dimensional lattice, known exact 
results for the asymptotic scaling of Z with the number of vertices can be expressed in a unified way 
as different values of one and the same function. 

1 Introduction 

Some problems in classical statistical mechanics have been found to have alternative formulations involv- 
ing fermionic (i.e. anticommuting) degrees of freedom. The most well-known examples are the Ising and 
dimer models [1] (which are related, as the former model can be mapped to the latter [2]). For planar 
graphs the generating function for the dimer model can be written as a Pfaffian [3], which can be regarded 
as a partition function for a system of noninteracting fermionic degrees of freedom [4]. The Pfaffian 
method can be generalized to nonplanar graphs of genus g, for which the dimer generating function can 
be written as a linear combination of 4 s Pfaffians, as first posited in Ref. [5], where this solution was 
explicitly demonstrated for the case of a square lattice embedded on a torus (g = 1); more general and 
rigorous discussions and proofs have later been given [6, 7, 8, 9, 10]. The Pfaffian method allows for 
a straightforward solution of the Ising and dimer models for two-dimensional lattices without crossing 
edges, as then the partition function involves only one or at most a few Pfaffians, and explicit expres- 
sions for these can be easily found using Fourier transformation provided the parameters describing the 
spin-spin interactions/dimer weights have some translational invariance. 

More recently, some other classical spin models [1 1], as well as generating functions for certain types 
of forests and trees (some of which are related to the q — > limit of antiferromagnetic Potts models) 
[12], have also been shown to have fermionic formulations. Notably, for these problems the fermions are 
interacting, so that even if one limits consideration to planar graphs, finding explicit exact solutions is not 
expected to be possible in general. 
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For all problems mentioned above, the fermionic formulations are most naturally given in terms of 
integrals over Grassmann variables [4] living on the vertices of some graph. A Pfaffian (whose square is 
a determinant) can be expressed as a Gaussian Grassmann integral, i.e. its "action" is a quadratic form in 
the Grassmann variables, corresponding to noninteracting fermionic degrees of freedom. In contrast, in 
the interacting case the action is non-Gaussian (typically quartic). 

In this paper we find fermionic formulations of another class of problems: the enumeration of g-edge- 
colourings of simple regular graphs of degree q. If each edge of a graph is coloured with one out of q 
possible colours, such that no edges connected to the same vertex have the same colour, the result is said 
to be a g-edge-colouring of the graph. A simple regular graph of degree q (from now on referred to as a 
graph of degree q for short) is a simple graph for which all vertices are connected to exactly q edges. We 
also note here that the ^-edge-colouring problem on a graph can be reformulated as a g-vertex-colouring 
problem (i.e. a zero-temperature g-state antiferromagnetic Potts model) on a closely related graph (see 
e.g. [13]). 

We first show that the number of g-edge-colourings Z of a graph of degree q can be obtained from 
the dimer generating function 3? for the same graph by successively differentiating 2£ q with respect to 
the dimer weight on each edge. This expression for Z can alternatively be rewritten in terms of connected 
dimer correlation functions. The Grassmann integral representation for Z is obtained by invoking the 
Pfaffian solution for H£ expressed in terms of Grassmann integrals. We consider planar graphs (including 
graphs that can be embedded on a cylinder) as well as (nonplanar) graphs that can be embedded on a 
torus. The Grassmann integral expressions obtained for Z for these two types of graphs differ due to the 
difference in the forms of the Pfaffian solution for 2£ ', the origin of which is the different topology of the 
embedded graphs, as mentioned above. 

Our final Grassmann integral expressions for Z (Eqs. (13) and (18)) involve q Grassmann variables 
(one for each colour) on each vertex and 2 Grassmann variables on each edge. Each term in the action 
couples the Grassmann variables on an edge to the Grassmann variables with a given colour index on the 
two vertices connected to the edge. Thus the action is spatially local and quartic. The action does not 
couple different colours directly; instead there is an indirect coupling via the edge Grassmann variables. 

The Grassmann integral formulation could serve as a starting point for further work based on analytical 
or numerical approaches. As an example of the latter, we discuss numerically exact evaluations of the 
Grassmann integrals using an algorithm by Creutz [14]. Applying this approach to the enumeration of 
4-edge-colourings on finite-size square lattices embedded on a torus, we verify that the results are in 
agreement with numerical transfer matrix calculations and also consistent with an earlier Bethe ansatz 
prediction for the thermodynamic limit [15]. Another possible approach (not pursued here) might be to 
make use of Hubbard-Stratonovich or related transformations (see e.g. [16]) to map the fermionic problem 
into a bosonic one (i.e. one involving integrals over ordinary c-numbers), which could then be analyzed 
using various methods. 

Exact results for the asymptotic exponential growth of Z with the number of vertices N have previously 
been derived for the honeycomb [17], square [15], 3-12 [13], and 4-8 [13] lattices. In an Appendix we 
show that the results for the first three of these, and for the one-dimensional chain, can be written in a 
unified way as different values of one and the same function. This suggests that a unified derivation of this 
asymptotic growth might be possible at least for this collection of lattices. We hope that the Grassmann 
integral formulations developed in this paper might be useful for shedding further light on this issue. 

This paper is organized as follows. Sec. 2 discusses the dimer formulations, with a technical derivation 
relegated to Appendix A. The fermionic (Grassmann integral) formulations are discussed in Sec. 3. The 
numerical evaluation of the number of 4-edge-colourings on the toroidal square lattice is presented in Sec. 
4, with some remarks on the computer implementation in Appendix B. Concluding remarks are given in 
Sec. 5. Appendix C presents the unifying formula for previous asymptotic results for Z for some lattices. 
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2 Dimer formulations 



2.1 Formulation in terms of the dimer generating function 

Consider an undirected graph G = (V,E) where V is the set of vertices and E the set of edges. We assume 
that the graph is simple (i.e. each edge connects two different vertices and any two vertices are connected 
by at most one edge) and that all vertices have the same degree q (i.e. each vertex has exactly q edges 
connected to it). Thus G is a simple regular graph of degree q; for short, we will refer to G as a graph of 
degree q. 

Labeling the edges as a = 1, . . . , \E\ where \E\ is the total number of edges, a dimer covering n of G 
can be represented as n = (n\ ,ri2, ■ ■ ■ where n a (= or 1) is the number of dimers on edge a, subject 
to the constraint that exactly one of the edges connected to each vertex hosts a dimer. Letting w a denote 
the weight of a dimer on edge a, the dimer partition function (generating function) 3f is 

n a 

Evaluating 3f with all dimer weights set equal to 1 gives the total number of dimer coverings. 

Next, we introduce coloured dimers and define a c-coloured dimer covering to consist of dimers that 
all have the same colour c, which may take one out of q values: c = l,2,...,q. Furthermore, we define 
a g-dimer-covering of G to be a composite structure of q dimer coverings, one of each colour. Because 
all vertices have degree q, any g-edge-colouring of G is also a g-dimer-covering, but unlike generic q- 
dimer-coverings it satisfies the additional "colouring constraint" that no edge should have more than one 
dimer, or, equivalently, no edge should lack a dimer. Thus the set of g-edge-colourings is a subset of the 
set of g-dimer-coverings. We therefore want to extract the former from the latter. To this end, consider 
the generating function for g-dimer-coverings, which is given by 

«('),...>) « 

where denotes a c-coloured dimer covering. Each term in the sum corresponds to a g-dimer-covering. 
A term that also corresponds to a ^-edge-colouring will contain exactly one factor of w a for each edge 
a. Conversely, for a term that does not correspond to a ^-edge-colouring, at least one edge a will contain 
more than one dimer and thus comes with a factor w p a with p > 2, and, equivalently, at least one edge j3 will 
not have a dimer and thus comes with a factor = 1 . The latter property implies that if we successively 
differentiate a term in (2) with respect to each dimer weight w r , the final result will be 1 for any term 
representing a ^-edge-colouring and for all other terms. Therefore the number of ^-edge-colourings Z 
is given by 




We have thus shown that the number of g-edge-colourings Z can be obtained from the generating function 
^ for the dimer problem with edge-dependent dimer weights on the same graph. 

2.2 Formulation in terms of connected dimer correlation functions 

In this subsection we will see that the number of ^-edge-colourings Z can be expressed in terms of con- 
nected correlation functions for the dimer problem. (Readers who are mainly interested in the Grassmann 
integral formulations of Z may skip this subsection, as later sections do not depend on it.) 
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The m-point dimer correlation function for the m edges a,\ , . . . , a m (where m can take values from 1 to 
\E\) is given by 



1 
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(4) 



The connected m-point dimer correlation function instead involves derivatives of log 3f in the usual way: 



,(m) 



(ai,...,<Xm) 



logiT. 



(5) 



(The subscript c on means 'connected' and should not be confused with the colour label introduced 
earlier.) These correlation functions are invariant under permutations of the m edge indices oti , . . . , a m . 
The correlation functions also depend on the dimer weights of all edges in the system, but (just as for the 
generating function £F) we do not indicate this dependence explicitly. 

Next note that the rhs of (3) is independent of the values of the dimer weights, so we can choose these 
at our convenience. By choosing to evaluate the rhs at w a = 1, we can replace d/dw a by w a d/dw a = 
d /d\ogw a . Thus Z can be rewritten as 



n - 

L a L dlogW a 



3f q 



(6) 



w a =l 



Starting from this expression, we show in Appendix A that Z can be written as 

B \E\ N P 

p : '- 1 ,„ I 



(7) 



The sum is over the B\ E \ partitions of E (here B? is the Bell number, the number of partitions of a set of 
I objects), and Np is the number of subsets in partition P. The integer r = l,...,Np labels the subsets 
Sp r of the partition, mp r is the number of edges in Sp r , and ^ mpr ' 1 (Sp r ) is the connected dimer correlation 
function for this subset. 



3 Fermionic (Grassmann integral) formulations 

In this section we derive fermionic formulations, in the form of expressions involving Grassmann inte- 
grals, for the number of g-edge-colourings Z of graphs of degree q that can be embedded without edge 
crossings either on a plane (Sec. 3.1) or on a torus (Sec. 3.2). It is instructive first to consider some exam- 
ples of such graphs, shown in Fig. 1. Graphs (a)-(d) are in the planar category while (e) is in the toroidal 
category. While the fermionic formulations are applicable to graphs of arbitrary size, including very small 
ones like graphs (a) and (b), regular lattice graphs of large (macroscopic) size are typically of greatest 
interest in statistical mechanics. Examples of such graphs in the toroidal category can be constructed very 
easily and naturally from Archimedean tilings by imposing periodic boundary conditions (BC's) in both 
directions. These tilings include many of the most frequently studied two-dimensional lattices in statisti- 
cal mechanics, including the square lattice on which graph (e) is based. Graphs (c)-(d) are two examples 
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of regular lattice graphs in the planar category that are most naturally embedded on a cylinder (rather than 
on a plane), which is also how they are drawn here (with open BC's in the horizontal direction and peri- 
odic BC's in the vertical direction). Graph (d) also exemplifies how a regular lattice graph in the planar 
category can be constructed by embedding an Archimedean tiling (in this case a honeycomb lattice) on a 
cylinder by imposing the appropriate BC's, and then introducing additional edges along the boundary (the 
vertical edges here) to ensure that the boundary vertices have the same degree as those in the bulk. 




(b) (c) (d) (e) 

Figure 1: Some examples of graphs of degree q (with q = 2,3,4) to which the methods developed in this 
paper can be applied to calculate the number of ^-edge-colourings Z (all graphs shown have Z / 0). In 
these figures, "loose" edges on opposite sides of each of the two directions (horizontal and vertical) should 
be identified. Graphs (a)-(d) are planar, with (c)-(d) drawn here in terms of their alternative (and more 
natural) embedding on a cylinder, while (e) is embedded on a torus. 



3.1 The planar case 

We start by considering a graph G = (E,V) of degree q as defined in Sec. 2.1. We may assume that the 
number of vertices N is even, as otherwise the possibility of dimer coverings of G (and in turn ^-edge- 
colourings) is trivially ruled out. We introduce an orientation of G defined by a set of arrows, one on each 
edge, each arrow pointing towards one of the two vertices connected to the edge. To a given orientation 
we associate a unique antisymmetric N xN matrix M with matrix element Mjj = if vertices i and j are 
not connected by an edge, and |M,y| = 1 if they are, with M ( y = +1 (—1) if the arrow on the edge points 
from i to j (from j to i). We will refer to M as a signed adjacency matrix. 

We next make a further restriction to planar graphs, which by definition can be embedded in the plane 
without crossing edges. From now on we consider such an embedding, which defines a set of faces. 
An orientation of the graph that has an odd number of edges oriented clockwise around each face will 
be called a Kasteleyn orientation (KO). Given a KO we can generate other KO's from it by Z2 gauge 
transformations, which conserve the clockwise-odd constraints. The elementary Z2 gauge transformation 
G, reverses the arrow directions on all edges connected to vertex i. A general Z2 gauge transformation is 
then given by G$ = HieS^i where S C V. Different KO's that are related by a Z2 gauge transformation 
are said to be gauge-equivalent. For planar graphs it turns out that all KO's are gauge-equivalent. Let 
K denote the signed adjacency matrix corresponding to a specific KO, and let A be the antisymmetric 
N x N matrix defined by A i; - = WijKij, where w,j (= Wp) is the edge-dependent dimer weight if i and j are 
connected by an edge, and zero otherwise. The dimer generating function is then given by [3] 

2T=\PfA\ (8) 

where Pf denotes the Pfaffian. As G, changes Pf A by a sign, Z2 gauge transformations leave 2? invariant. 
If all dimer weights are set to unity, A reduces to K and 3? becomes the number of dimer coverings. By 
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introducing a Grassmann variable l//,- on each vertex i, Pf A can be expressed as a Gaussian Grassmann 
integral over these variables: 1 



Pf A = J d\\r N ■ ■ -dyi exp ^ £ Yi A ijYjj ■ 



Thus 



3f q = 



I 



(c) 



c=l 



ex p( ^LLW'W^j 



c=l ij 



(9) 



(10) 



where we introduced a colour superscript c = 1 , . . . , q to distinguish different Grassmann variables be- 
longing to the same vertex. Using Eq. (3) then gives 
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(11) 



(In the product over edges, we have labeled the two vertices touching an edge a as i\ a and i2 a . Thus in 

(c) 

this product, a given Grassmann variable x//- will in fact appear under q different names, one for each 
edge the vertex i touches.) This expression can be simplified further. First, note that expanding out the 
product over edges gives a sum of terms, each of which contains 2\E\ = qN Grassmann variables (not 
necessarily all distinct), which also equals the total number of Grassmann variables being integrated over. 
From the properties of Grassmann variables it then follows that the exponential can be replaced by 1 . This 
shows that the rhs of Eq. (11) is independent of the dimer weights w a . Second, the signs Ki la j 2a in the 
product over edges are superfluous and can be omitted. Thus 



/ 



(c) 



n y v- (Ca V Co) 
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a c„=l 



(12) 



Next we introduce a pair of Grassmann variables E, a and E, a on each edge a (for an edge connecting 
vertices i and j we can write the variables as = and fyj = ^ 7I ). Then Z can be written as 



Z = 



(c) 



c=\ 



UdSadSa exp I ^IlWX4^^j 



■ c=\ ij 



(13) 



Here M is the signed adjacency matrix corresponding to an arbitrary orientation of the graph. To prove 
this result one simply integrates out the edge Grassmann variables, which leads back to Eq. (12). 

Eqs. (12) and (13) are the main results of this subsection. Although Eq. (13) does in some sense 
represent a less economical formulation than Eq. (12) due to the additional Grassmann variables living on 
the edges, Eq. (13) is of interest because it takes the "standard" form /[integration measure] exp(action) 
for a partition function, with the action a sum of spatially local terms; this form may be more convenient 
for further manipulations. The fact that the action in Eq. (13) is quartic in the Grassmann variables, i.e. 
non-Gaussian, is a reflection of the interacting nature of the problem. We also emphasize that since M in 
Eq. (13) can be any signed adjacency matrix for the graph, Eq. (13) for Z is invariant under reversals of 
the orientations of an arbitrary subset of edges. This is a much larger set of transformations than the Z2 
gauge transformations that leave the dimer generating function 5° invariant. 

The validity of Eq. (13) can be verified using a somewhat different line of reasoning, starting from the 
observation (cf. the discussion in Sec. 2.1) that the calculation of the number of ^-edge-colourings Z is 



1 An equivalent representation is Pf A = / d . . . d Xj/^ exp ( — \ Wi^ij Wj ) 
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similar to the calculation of the number of g-dimer-coverings 2£ q \ Wa= \ , but with the crucial difference that 
all g-dimer-coverings that contain edges with multiple dimers are not included in the count. This suggests 
that in order to find Z one can first modify Eq. (10) (with all dimer weights set to unity) by making 
the replacement Kjj — > Kij = Kij^ij^j. We can think of the Kjj as matrix elements of a "Grassmann- 
valued Kasteleyn matrix". As these matrix elements contain a product of an even number of Grassmann 
variables, they behave like ordinary c-numbers (in particular, they commute with everything) except that 
their square vanishes by virtue of f = = 0. This latter property ensures that all g-dimer-coverings 
containing edges with multiple dimers will give zero contribution. Integrating over the vertex Grassmann 
variables therefore gives \ZY\ a E, a E, a \. Introducing also an integral over the edge Grassmann variables, 
one is left simply with Z. Finally, one can convince oneself that the same final result is obtained if K is 
replaced by an arbitrary signed adjacency matrix M for the graph. 



3.2 The toroidal case 



Here we consider nonplanar graphs of degree q with the same defining properties as those discussed in 
Sec. 3.1, except that these graphs can be embedded on a torus, not on a plane (without crossing edges). As 
before, an orientation is described by a signed adjacency matrix M. Kasteleyn orientations and Z2 gauge 
transformations are also defined in the same way as for planar graphs. For toroidal graphs the Kasteleyn 
orientations fall into 4 gauge-inequivalent classes. The dimer generating function 3? can be written 



£ r M Pf A<") 
li=l 



(14) 



where the sum goes over the 4 gauge-inequivalent classes, A^ = w^K^, where is the Kasteleyn 
matrix corresponding to a Kasteleyn orientation chosen from class jl, and the = ±1 are appropriately 
chosen signs. We refer to the literature for more detailed discussions of these matrices and signs [5, 6, 7, 
8,9, 10]. 

Following the same steps as in Sec. 3.1, we now find 



and 
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Z = 2~ q 



Ml Pel 



I 



Hdy^ ...d\j/\ c 



c=l 



■ c=l ij 



(15) 



(c) 



c=l 



exp 
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a c a =l 



a, '2a 



(16) 



The exponential can be replaced by 1 by the same argument as for the planar case. The signs contributed 
by the signed adjacency matrices can however not be completely eliminated, unlike the situation for the 
planar case. But Z is still invariant under reversal of the orientations of an arbitrary subset of edges. As 
noted before, this set of transformations is much larger than the Z2 gauge transformations, and thus the 
signed adjacency matrices resulting from these transformations will generally not satisfy clockwise- 
odd constraints around the faces. Note however that the change in is independent of jl and thus the 
transformations preserve the ratios Afjf /Atjf\ We get 



Z = 2 



I 



.dYi 



n 1 <>, 



(ca)*Af*c a ) 
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ha,i2a Vha 



(17) 
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Note that Af,^,- is independent of }i for some subset of the edges. For these edges this sign is superfluous 
in this expression and may be omitted. 

Next, by introducing Grassmann variables f a and E, a on each edge a, we can express Z as 



Z = 2~ q 



■r. 



(c) 



=1 «/ 



(18) 



Eqs. (17) and (18) are the main results of this subsection. They are the torus analogues of Eqs. (12) and 
(13) for the planar case. When properly generalized, the remarks made at the end of Sec. 3.1 are valid also 
for Eqs. (17) and (18). The only difference between the two cases is that in the toroidal case the analysis 
involves four different signed adjacency matrices instead of one, which leads to a sum of 4 q terms 
in the expression for Z. The number of distinct terms is however less than 4 q due to the invariance of the 
summands in Eqs. (17) and (18) under permutations of fX\, . . . ,n q . Moreover, as we will see an example 
of in the next section, additional simplifications may be possible which can allow one to further reduce, 
to a handful or less, the number of terms that need to be explicitly evaluated. 



4 Exact numerical evaluation of Grassmann integral expressions: appli- 
cation to 4-edge colourings on a square lattice embedded on a torus 

In this section we present an approach for explicit calculation of the number of ^-edge-colourings Z on 
graphs of degree q, based on numerically evaluating expressions for Z like Eqs. (13) or (18). These involve 
Grassmann integrals on the form / dr\ n . . . Jrji exp(5({Tj })) which can be evaluated using an algorithm due 
to Creutz [14]. 

As a concrete example of this approach, we consider the 4-edge-colouring problem on a toroidal 
square lattice (cf. Fig. 1(e)) with N x (N y ) vertices in the x (y) direction. Thus Z can be found from Eq. 
(18), with appropriate choices for the matrices AfM and the signs obtained from the Pfaffian solution of 
the dimer problem on this lattice [5]. While Eq. (18) in this case involves a sum of 4 4 = 256 Grassmann 
integrals, the invariance of the summand under permutations of jUi , - - - , /X4 can be used to reduce the 
number of terms to 35. 

Some aspects of our implementation of Creutz's algorithm are discussed in Appendix B. The nature of 
the algorithm and the specific problem to which it is applied, together with generic limitations on computer 
speed and memory, imply that in practice the maximum value of N x that can be considered is rather small; 
in our implementation we were able to go up to N x = 8. The maximum values of N y can be considerably 
larger but decrease as a function of N x . For all calculations done with N x = 3-6 we evaluated all 35 terms 
and observed that many of them evaluated to zero while many others cancelled among themselves. As a 
consequence the correct answer for Z could in fact be obtained by restricting evaluation to a much smaller 
number of terms: 5 terms if both N x and N y are even, and only 1 term if either N x or N y is odd. This 
made it possible for us to also calculate Z for cases with N x = l and 8 that otherwise would have been 
too time-consuming. Table 1 shows our results for Z = Z(N x ,N y ) for N x = 3-8. Note that for N x = 8, we 
restricted the largest N y values to be odd, as this allowed us to find Z by calculating only one Grassmann 
integral. 

To check the correctness of the results in Table 1 we have compared them in various ways to results 
from numerical transfer matrix calculations and also to a previous Bethe ansatz study [15]. To this end, 
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UO OU7 1UO Z.U J U7U OOZ. UUJ UZ. 1 / 7Z / UO 


4x10 


67 658 756 896 


7x8 


1 1 607 1 47 777 448 

1 1 UU / It / Z. / Z. fro 


4x11 


742 987 407 360 


7 x 10 


1 486 277 1 24 7 1 5 397 

±U*TOUZ./ / IZt / 1J J7i 


4x12 


8 Q4? 459 8QQ 840 
o ?4z hjz oyy otu 


7 x 12 


10 ?17 845 ?Q1 845 1 Q5 808 

lUZJ/ Ot J z.7 1 o^t J 1 7J ouo 


4x 13 


106 992 869 867 520 


7 x 14 


10 322 131 976601 474 883 392 


4x14 


1 7X4 867 476 349 568 

1 Z.OT OUZ, T-Z.U J4 " JUO 


7 x 16 


1 557 848 759 1 99 993 065 837 056 

1 U J JZ OHO 177 77J UUJ O J / U JU 


4x 15 


15 407016206008 320 


7 x 18 


10855518364 633 186 870 098 974592 


4 x 16 


184918 168859 836416 

1 ~J 1 O 1UO O J 17 O J U 1 U 


7 x 20 


11 1 98 03 1 662 235 847 304 452 806 203 520 

1 ± 170 U J 1 OUL iJJ O^T / JU^T T J L OUU ZU J ~J £J\J 


4x17 


? ?1 8 61 1 090 841 61 5 160 

AAlOUl 1 UZ.U Oi 1 Ul J JUU 


7 x 22 


1 1 566 056 69Q 6Q1 191 689 661 1 95 ?45 140 41 6 

1 1 JUU U JU UZ7 u7l 1Z1 UOZ. UUJ 1ZJ i^J Jt-U TlU 


4x18 


26 624 552 900 1 85 030 656 

j-. u ult jjl y\/\j ioj u ju u j u 


7 x 24 


1 1 953 1 67 548 050 678 272 222 497 989 683 044 256 

1± 7J J 1U/ JtOUJUU/OZ/LLZZ*t;// 7u / OOJ Ut7 LJU 


4x19 


11 Q47Q QQ7 QQ6 915 48Q 980 

J 17 *T / " "y 1 77U liJj tO? iOU 


7 x 26 


1 ? 156 574 8Q1 854 1 11 6?6 ?Q7 HQ 854 605 QO? 1 ?5 948 

JJU J/TU71 OJ4 1 JJ UZ.U Z." / JJ7 OJl UUJ 7UZ. 1 Z. J iiu 


4 x 20 


1 811 801 880 6Q4 041 1 1 5 590 

J O J J OUJ OOU Kjy^r UT'J 1 1 J JZU 


7 x 28 


1 9 775 1 Q7 Q10 055 919 189 651 748 518 401 498 081 77Q 198 

1Z / / J lyl 7JU UJJ Z1Z JOZ. U J 1 / T^O J JO tuJ tZO UOl / / 7 JZ.O 


5x6 


1 7 Q88 Q60 

i / yoo 7UU 


7 x 30 


1 1 908 779 604 91 1 Q70 774 HQ 1 Ql 141 67Q 044 41 9 01? 615 1 68 
i j z.uo / / z. uut z.ii y i\j i /t JJ7 i"J iii \j i y Tii ujz ujj i uo 


5x8 


1 840646400 


7x32 


13 657 432 485 307 219 457 227 579 236 006 517 304 339 033 098 752 


5 x 10 


209 977 817280 


8x8 


1685 928 423 086592 


5 x 12 


24 836 803 964640 


8x9 


65 640522173 916672 


5 x 14 


2 974 334 794053 120 


8x 11 


361943 772 232 276 810752 


5 x 16 


357 739836702 854400 


8x 13 


1 983 974 1 00 428 34 1 796 935 680 


5 x 18 


43 094084170133 825760 


8x 15 


10 860 078 679 370 7 14 987 464 273 920 


5x20 


5194113814423 956157 440 


8x 17 


59 428 139 978 826 027 050 009 652 486 144 


5x22 


626 172 572 098 389 717 579 840 


8x 19 


325 176 895 253 864 837 579 133 026 1 12 749 568 



Table 1: The number of 4-edge-colourings Z = Z{N x ,N y ) on a square lattice of size N x x N y embedded on 
a torus. The results were obtained by exact numerical evaluation of Eq. (18), using the algorithm in [14] to 
evaluate the Grassmann integrals. Only cases with < are shown, as Z is invariant under interchange 
of N x and N y . (For N X = S with N y > 9 only odd values of N y were considered, as the computation times 
for even values of N y were prohibitively large.) 

we write Z as 

Z = Tr T N > ' = £Af- (19) 

where the sum is over the eigenvectors i of the 4 Nx x 4 Nx row transfer matrix T with eigenvalues A;. We 
have diagonalized T numerically for N x = 3-6 and verified that evaluating Z from Eq. (19) reproduces the 
results in Table 1 for not too large values of (this maximum value of decreases with A*; for larger 
values of Ay computational errors due to the floating-point arithmetic involved in the transfer matrix 
calculation of Z become too large). Our results for A* = 7, 8 with A^ < 6 are also indirectly verified in this 
way due to the symmetry Z(N x ,N y ) = Z(N y ,N x ). 

For our further checks it is useful to consider f(N x ,N y ) = (N x N y )~ l logZ(N x ,N y ) which in Fig. 2(a) 
has been plotted as a function of N~ l for fixed values of A* (with N y restricted to values such that 
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N x 


^max 




f(N x ,oo) 


^max 


3 


4 


4 


0.41415 


3.4641 


4 


1 





0.62123 


12.000 


5 


4 


4 


0.47922 


10.980 


6 


1 





0.55919 


28.649 


7 


4 


4 


0.49580 


32.156 


8 


1 





0.53796 


73.971 



Table 2: Values of <i max , d min , f{N x ,°°), and = exp(N x f(N x ,°°)) obtained by fitting f(N x ,N y ) in Fig. 
2(a) to the expression for f(N x ,N y ) predicted by the transfer matrix approach for large N y . For N x = 3- 
6 we have verified the values listed here by also finding A max , <i max , and <i m i n directly from numerical 
diagonalization of the transfer matrix. 

N x N y is even, as Z = otherwise). The oscillating behaviour of / for even N x can be understood as 
a consequence of negative eigenvalues of T, as these contribute to Z with opposite signs for even and 
odd N y . As N~ becomes small, / is seen to approach a straight line in all cases. To understand this 
we first note that for large N y the sum in Eq. (19) is dominated by the terms involving the eigen- 
values with largest modulus, i.e. Z ~ A™ x [fi? max + d m { n (— l) Ny ], where Am ax > is the largest eigen- 
value, <i max is its multiplicity, and d m i n is the multiplicity of the eigenvalue — Am ax . Thus in this limit 
f(N x ,N y ) ps N~ l log Am ax +N~ l N~ 1 log[<i max + (-l) Ny d m i n ]. By fitting this expression to the small-A 7 " 1 
region of the curves in Fig. 2(a) one can deduce the values of f(N x ,°°), Am ax , d max , and d mm . The results 
are given in Table 2. Our numerical diagonalization of T for N x = 3-6 reproduces these values for those 
cases. 

We briefly discuss how the fitting is done. First consider the case of odd N x , for which the slope of 
the straight line is N~ 1 log [<i max + <i mm ]. Since d max + d mm must be an integer, its value can be deduced 
with certainty from the slope, and in turn the exact value of the slope can be found. Next one can estimate 
the intercept f(N x ,°°) with the vertical axis, from which one gets A max = exp(N x f(N x ,°°)). Finally, since 
Z = for odd N y it follows from the large-A^ expression for Z in this case that <i m i n — <^max- For the 
case of even N x the deduction of the multiplicities is slightly different. From Fig. 2 one sees that in this 
case / becomes increasingly independent of N~ 1 as N~ 1 is decreased. Comparing this behaviour to the 
coefficient N~ l log[d max + (—l ) Ny d m i n ] of the Ny 1 term in /, one is able to conclude that d max = 1 and 
dmin = (and so the slope is identically zero). 

Finally, we show in Fig. 2(b) f(N x ,°°) as a function of A^T 1 together with the Bethe ansatz prediction 
[15] for /(°°,°°) (horizontal dashed line), given by log(r 2 (l/4)/[V2^ 3 / 2 ]) « 0.51238. The observed 
behaviour of f(N x ,°°) as A^" 1 decreases is clearly consistent with this prediction. 

5 Concluding remarks 

We have shown that the number Z of ^-edge-colourings of simple regular graphs of degree q are deducible 
from the dimer generating function 3f (or, equivalently, from the set of connected dimer correlation 
functions) on the same graph. For graphs of this type that are either planar or embeddable on a torus we 
invoked the expressions for 3f in terms of Pfaffians of Kasteleyn matrices to derive fermionic expressions 
for Z in the form of Grassmann integrals. In particular, we derived expressions (Eqs. (13) and (18)) in 
which the integrands are given as the exponential of a non-Gaussian quartic "action" that is a sum of 
spatially local terms involving both vertex and edge Grassmann variables. 
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Figure 2: (a) f(N x ,N y )= (N x N y ) 1 logZ(N x ,N y ) plotted as a function of N y 1 for fixed values of N x , with 
Z obtained from Table 1. (b) f(N x ,°°) as a function of A^T 1 (symbols are the same as in (a)). The values 
of f(N x ,°°) (listed in Table 2) were obtained as the intercepts of the extrapolation of f(N x ,N y ) in (a) with 
the vertical axis N~ l = 0. As A^T 1 — > we expect f(N x ,o°) to approach the Bethe ansatz prediction [15] 
for /(°°, o°) (the horizontal dashed line). In both (a) and (b), the straight lines connecting points are guides 
to the eye. (Such lines were not drawn for the N x = 8 case when N y > 9 as in this parameter regime Z was 
not calculated for even values of N y ; cf. remark in caption of Table 1.) 



As an example of a possible use of these expressions, we discussed exact numerical evaluations of 
them using an algorithm by Creutz [14], and presented a concrete application to the enumeration of 4- 
edge-colourings on a square lattice embedded on a torus. The agreement found with results from Bethe 
ansatz [15] and numerical transfer matrix calculations serves as a direct nontrivial check of the correctness 
of Eq. (18). We note that, unlike the numerical transfer matrix method, the numerical Grassmann integral 
approach only involves integer arithmetic, and thus numerical floating-point errors are never an issue. 
Another difference is that the Grassmann integral approach can also be applied to graphs without any 
repeated structure. (See however [18] for a discussion of a generalized transfer matrix approach that also 
has these properties.) 

Acknowledgements 

We acknowledge use of the C++ linear algebra library Armadillo [19] for the numerical transfer matrix 
diagonalizations discussed in Sec. 4. 

A Proof of Eq. (7) 

We start by defining the auxiliary quantity 

where £ = 1,2, . . . , \E\. Thus from Eq. (6) we have 

Z = ©(|£|)L=r (2D 
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Explicit evaluation gives 



0(1) = 


2? q i 


7^ (1) («i), 


0(2) = 


2? q 


V^c (1) («l 


0(3) = 


geq 


v^ (1) («l 



(22) 
(23) 



?(!)/ 



?(!)/ 



(«! , a 2 )^ a) (a 3 ) + ^' («i , a 3 )9T (« 2 ) 



(2), 



+ ^ 2) (a 2 ,a 3 )^ lj (ai))+^^(ai,a 2 ,«3) . 
From these expressions a pattern can be discerned, which suggests the following ansatz for &(£): 

B t N P t 

i ^ n 

P,=\ r,=l 



7(3)/ 



(24) 



©(£)=^£^n^ (mW ( 5 ^)- 



(25) 



The sum is over all partitions of the edge set {(X\ , a 2 , . . . , a?} = which we label P( = \ ,2,. . . ,B(, where 
the Bell number Bf is the total number of partitions. Furthermore, ri = 1,2, . . . labels the subsets Sp m 

of partition Pi, with the total number of subsets, mp t . rf the number of edges in Sp er( , and ( $ < f lPlH ' \Sp in ) 
the connected dimer correlation function for this subset. 

We will prove (25) by induction. Assuming it holds for £ edges (the set E(), we consider the case of 
I + 1 edges (the set Ei + \ =EgU {(X£ + \ }). We have 



&(£+\) = 



d log w, 



-&(£) = 



<xe+i 



dlogw am 

N Pf 



r> I q Np < n ■ ' 

p t =l r t =\ 



(m Pf r f ) 



Pen 



(5p^) 
5 



z 

P/=1 



— n< 



(Sp tn ) 



X* £ ^ rfWi) U^ mpert \sp m ) 



r e =l 



Np f N P( 

+ £ ^ +1) (%< u{a, +1 }) n % (mPn \sp tn ) 



(26) 



This expression is in fact of the form (25) with £ replaced by £ + 1 . To see this, first note that the partitions 
of E( + i can be generated by starting from those of Eg and "adding" the additional edge in all possible 
ways. For a given partition Pi of Ei we can either let the additional edge be a subset by itself [this 
corresponds to letting <9/dlogw a/+1 act on 3f q , which gives a factor 3? q (o^+i); this produces the 
first term in Eq. (26)], or we can add it to any one of the Np ( subsets in that partition [this corresponds 
to letting d /d\ogw ae+l act on the connected dimer correlation function for that subset, which gives the 
connected dimer correlation function for the union of the original subset and the new edge; this produces 
the second term in Eq. (26)]. It follows that the total number of partitions Bi + \ of the set Ef + i is given by 



B t+1 = £(1+JVt%). 



(27) 



ft=i 



Eq. (27) can be verified by invoking the so-called Stirling number of the second kind, S{£,k), defined as 
the number of ways to partition a set of £ objects into k (non-empty) subsets. From this definition it follows 
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that the Bell number Be can be written B( = E|=o^(A&) an( ^ furthermore that Yfp t= \Np t = Y,l = okS(£,k). 
Eq. (27) can then be shown to follow from the recurrence relation S(£ + l,k) = S(£,k — 1) + kS(£,k). 

Thus we have proved that if Eq. (25) holds for £, it also holds for £+ 1. Since Eq. (25) is true for 
£=l, the induction proof of Eq. (25) is complete. Eq. (7) then follows from Eq. (21). 

B On our implementation of Creutz's algorithm 

In this Appendix we make some remarks on our implementation of Creutz's algorithm for the numerical 
evaluation of the Grassmann integrals in Eq. (18) for the problem of a square lattice embedded on a torus 
discussed in Sec. 4 (these remarks are not likely to be very comprehensible unless the reader is already 
familiar with the discussion of the algorithm in Ref. [14]). 

The x (y) direction is taken to be the transverse (longitudinal) direction in the algorithm. Memory 
requirements grow exponentially with N x , the number of vertices in the transverse direction. Moreover, 
the growth rate for our problem is large due to (i) the large effective number of Grassmann variables per 
vertex, and (ii) the periodic boundary conditions in the longitudinal direction (which effectively doubles 
the growth rate compared to the case of open boundary conditions [14]). Efficient use of memory is 
therefore important. We now describe some simple ways to reduce the memory needed. 

Regarding (i): The effective number of Grassmann variables per vertex is naively 8, since there are 4 
Grassmann variables living on each vertex, and 2 Grassmann variables living on each of the 2 edges that 
can be associated with a given vertex. The algorithm introduces one fermionic single -particle state for each 
Grassmann variable. Applying the algorithm to our problem it can however be seen that the fermions in 
the two single-particle states on a given edge are always annihilated together. It is then possible to modify 
the algorithm slightly so that for each edge one only needs to introduce a single fermionic single-particle 
state. This reduces the effective number of Grassmann variables per vertex to 6 (which is still large). 

Regarding (ii): The memory requirements can be considerably reduced by splitting the calculation of 
the Grassmann integral up into different parts. Let v = 1 , . . . , 6N refer to a Grassmann variable (N = N x N y ). 
The variables are labeled according to the vertex they are associated with, with vertices 1 and N in the 
bottom left and top right corners respectively. The Grassmann integral can be written [14] (0| Y6n) where 
|0) is the state with no fermions and \\j/ v ) = <2 v |Vv-i), where Q v is an operator and \\j/o) = \F), the 
completely occupied state. Consider the state \Y6N X } obtained after integrating out all 6N X Grassmann 
variables associated with vertices in the bottom row. Each basis state with a nonzero coefficient in | Y6N X ) 
will have the following property: for each vertex in the top row (which is, due to the periodic boundary 
conditions, coupled to the vertex in the bottom row directly "above" it), one of the 4 fermionic states (each 
of which is associated with a different colour) will have occupation number 0, while the 3 other states will 
have occupation number 1. Thus one can split the terms in | Y6N X ) mto distinct parts characterized by which 
of these states have occupation number 0, i.e. one writes \Y6N x ) = Ep I Wf>N x ,p) where p labels the different 
parts. (The splitting can be done to various degrees, depending on one's needs: the coarsest splitting is 
into only 4 parts based on the states on just a single vertex in the top row, while the finest splitting would 
be into 4 Nx parts based on the states on all vertices.) Next one considers each part p in turn, calculating 
\Y6N,p) = (llv=6Af +1 Qv) I W(>N x ,p) where the £2v's are applied in ascending order of v. The final result is 
obtained as = £ P (0| WdN,p)- Doing this splitting has the advantage that the calculation of \ Y6N,p) 

involves a considerably smaller maximum size of the hash table than would a direct calculation of | Y6n) 
(i.e. without doing this splitting). This is because the calculations of | Y6N,p) for different p's have no 
states in common during those stages of the calculations when the hash tables reach their maximum sizes 
(which happens for specific intermediate values of v, not too close to the bottom or top rows). 

Each key-value pair in the hash table ("dictionary") consists of a many-particle fermionic basis state 
and its coefficient, both of which are integers (for the basis state, its binary representation gives the oc- 
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cupation numbers in each single-particle fermion state) whose size will typically exceed the computer's 
limit. Thus one needs a way to do arbitrary-precision integer arithmetic. We wrote the code in Python, 
which has built-in support for both hash tables and arbitrary-precision integer arithmetic. The code was 
run on various computers with 2.4-3 GHz processors and 4-16 GB RAM. For N x = 3 the calculation of a 
Grassmann integral in Eq. (18) took a few seconds or less, while for N x = 7,8 it took about 10 days for 
the largest values of N y considered. 



C Unified expression for W = Yim^^^Z 1 ^ for some lattices 

The problem of enumerating the g-edge-colourings for a lattice with degree (i.e. coordination number) q 
has been solved exactly, in the thermodynamic limit, only for a relatively small number of lattices: the 
honeycomb lattice [17], the square lattice [15] (see also [20]), the 4-8 lattice [13], and the 3-12 lattice [13]. 
These lattices have q = 3 except the square lattice which has q = 4. In addition the problem can be solved 
trivially for a one-dimensional lattice (q = 2). For a given lattice the exponential scaling of the number 
of ^-colourings Z is given by the parameter W = lim^^ooZ 1 /^, where N is the number of lattice sites 
(vertices). In the following we demonstrate a connection between the results for W for the honeycomb, 
square, 3-12, and the one-dimensional lattice. 

For the honeycomb (he) lattice, it was found that [17, 21] 

2 (3j-l) 2 _ 3r 3 (l/3) 

Whc "M^7^)"^^' (28) 

while for the square (sq) lattice, it was found that [15] 

- (4j-l)(2j-l) r 2 (l/4) 

Wsq = H 2j(4j-3) =VT^- (29) 

1 /3 

For the 3-12 lattice, W3-12 = W hc [13]. Finally, as the one-dimensional (ID) lattice has Z = 2 for all 
(even) N, it follows that Wi D = 1. 

The results for the parameter W for these four lattices can be written in a unified form in the following 
equivalent ways: 

W P = Y^{qj-l)(qj-q + 2) T[\/q) 



jJi qj(qj-q+l) T(l - l/q)T(2/q) 

sin(n/q) 



% 



'-B{\/q,\/q)= 2^(1/4,1 -2/4; 1;1), (30) 



where r(x) is the gamma function, B(x,y) is the beta function, and 2F1 (a, b;c;z) is Gauss's hypergeometric 
function. 2 Each lattice is associated with a pair of numbers (q,p) where q is the degree/coordination 
number as before and the value of p coincides with the number of sites per unit cell of the lattice. Thus 
(q,p)hc = (3,2), (q,p) sq = (4,1), (q,p) 3 _ u = (3,6), and (q,p)m = (2,1). 
On the other hand, for the 4-8 lattice, it was found that [13] 

00 / 1 _ v 4 ;'- 1 \ ! / 2 

2 In order to convert between the different expressions in (30) it is convenient to insert the factor T(l) = 1 in the numerator in 
the expression involving gamma functions. 
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where x = 2 — (p with <p = (1 + V / 5)/2 being the golden ratio. More recently this value for W4-& was 
also verified numerically [22]. However, using (q,p)4-g = (3,4) in (30) fails to reproduce this value. Of 
course, one might question whether the interpretation of p as the number of sites per unit cell is really the 
correct one. But keeping q = 3, numerical agreement between (30) and (31) would require one to take 
p « 1.759, a value that it is difficult to make sense of, also in light of the "natural" values for p required 
for the other lattices. Thus we conclude that (30) is not valid for the 4-8 lattice. 

A deeper understanding of (30) will have to await a general derivation of it for a class of lattices that 
should include the honeycomb, square, 3-12, and one-dimensional lattice, and possibly other lattices as 
well that have not been identified yet. 
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